Simulacion del gas ideal

Fuente: https://github.com/rafael-fuente/Ideal-Gas-Simulation-To-Verify-Maxwell-Boltzmann-distribution

Código de la simulación (tarda como 10 minutos)

In [ ]:
import numpy as np


class Particle:
    """Define physics of elastic collision."""
    
    def __init__(self, mass, radius, position, velocity):
        """Initialize a Particle object
        
        mass the mass of particle
        radius the radius of particle
        position the position vector of particle
        velocity the velocity vector of particle
        """
        self.mass = mass
        self.radius = radius
        
        # last position and velocity
        self.position = np.array(position)
        self.velocity = np.array(velocity)
        
        # all position and velocities recorded during the simulation
        self.solpos = [np.copy(self.position)]
        self.solvel = [np.copy(self.velocity)]
        self.solvel_mag = [np.linalg.norm(np.copy(self.velocity))]
        
    def compute_step(self, step):
        """Compute position of next step."""
        self.position += step * self.velocity
        self.solpos.append(np.copy(self.position)) 
        self.solvel.append(np.copy(self.velocity)) 
        self.solvel_mag.append(np.linalg.norm(np.copy(self.velocity))) 
        

    
    
    def check_coll(self, particle):
        """Check if there is a collision with another particle."""
        
        r1, r2 = self.radius, particle.radius
        x1, x2 = self.position, particle.position
        di = x2-x1
        norm = np.linalg.norm(di)
        if norm-(r1+r2)*1.4 < 0:
            return True
        else:
            return False

    
    def compute_coll(self, particle, step):
        """Compute velocity after collision with another particle."""
        m1, m2 = self.mass, particle.mass
        r1, r2 = self.radius, particle.radius
        v1, v2 = self.velocity, particle.velocity
        x1, x2 = self.position, particle.position
        di = x2-x1
        norm = np.linalg.norm(di)
        if norm-(r1+r2)*1.1 < step*abs(np.dot(v1-v2, di))/norm:
            self.velocity = v1 - 2. * m2/(m1+m2) * np.dot(v1-v2, di) / (np.linalg.norm(di)**2.) * di
            particle.velocity = v2 - 2. * m1/(m2+m1) * np.dot(v2-v1, (-di)) / (np.linalg.norm(di)**2.) * (-di)
            

    def compute_refl(self, step, size):
        """Compute velocity after hitting an edge.
        step the computation step
        size the medium size
        """
        r, v, x = self.radius, self.velocity, self.position
        projx = step*abs(np.dot(v,np.array([1.,0.])))
        projy = step*abs(np.dot(v,np.array([0.,1.])))
        if abs(x[0])-r < projx or abs(size-x[0])-r < projx:
            self.velocity[0] *= -1
        if abs(x[1])-r < projy or abs(size-x[1])-r < projy:
            self.velocity[1] *= -1.


def solve_step(particle_list, step, size):
    """Solve a step for every particle."""
    
    # Detect edge-hitting and collision of every particle
    for i in range(len(particle_list)):
        particle_list[i].compute_refl(step,size)
        for j in range(i+1,len(particle_list)):
                particle_list[i].compute_coll(particle_list[j],step)    

                
    # Compute position of every particle  
    for particle in particle_list:
        particle.compute_step(step)




################################################################################################################################





def init_list_random(N, radius, mass, boxsize):
    """Generate N Particle objects in a random way in a list."""
    particle_list = []

    for i in range(N):
        
#         v_mag = np.random.rand(1)*6 # CAMBIE ESTO
        v_mag = 3
#         v_ang = np.random.rand(1)*2*np.pi # CAMBIE ESTO
        v_ang = np.random.rand(1)*0.5*np.pi + np.pi
        v = np.append(v_mag*np.cos(v_ang), v_mag*np.sin(v_ang))
        
        collision = True
        while(collision == True):
            
            collision = False
#             pos = radius + np.random.rand(2)*(boxsize-2*radius) 
            pos = radius + np.random.rand(2)*(0.5*boxsize-2/2*radius) ### CAMBIE ESTA
            newparticle = Particle(mass, radius, pos, v)
            for j in range(len(particle_list)):

                collision = newparticle.check_coll( particle_list[j] )

                if collision == True:
                    break

        particle_list.append(newparticle)
    return particle_list



particle_number = 300
boxsize = 200.

# You need a larger tfin and stepnumber to get the equilibrium state. But the computation takes more time.
tfin = 130#*30
stepnumber = 2000#*30

timestep = tfin/stepnumber

particle_list = init_list_random(particle_number, radius = 1.5, mass = 1, boxsize = 200)

# Compute simulation (It takes some time if stepnumber and particle_number are large)
for i in range(stepnumber):
    print(i)
    solve_step(particle_list, timestep, boxsize)
    #print(i)


################################################################################################################################


# Visualization of the solution with matplotlib. It use a slider to change the time


# import matplotlib.pyplot as plt


# fig = plt.figure(figsize=(12, 6))
# ax = fig.add_subplot(1,2,1)

# hist = fig.add_subplot(1,2,2)

# plt.subplots_adjust(bottom=0.2,left=0.15)

# ax.axis('equal')
# ax.axis([-1, 30, -1, 30])

# ax.xaxis.set_visible(False)
# ax.yaxis.set_visible(False)
# ax.set_xlim([0,boxsize])
# ax.set_ylim([0,boxsize])

# # Draw Particles as circles
# circle = [None]*particle_number
# for i in range(particle_number):
#     circle[i] = plt.Circle((particle_list[i].solpos[0][0], particle_list[i].solpos[0][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
#     ax.add_patch(circle[i])

# # Graph Particles speed histogram
# vel_mod = [particle_list[i].solvel_mag[0] for i in range(len(particle_list))]
# hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
# hist.set_xlabel("Speed")
# hist.set_ylabel("Frecuency Density")



# Compute 2d Boltzmann distribution

#total energy should be constant for any time index
def total_Energy(particle_list, index): 
    return sum([particle_list[i].mass / 2. * particle_list[i].solvel_mag[index]**2  for i in range(len(particle_list))])


E = total_Energy(particle_list, 0)
Average_E = E/len(particle_list) 
k = 1.38064852e-23
T = 2*Average_E/(2*k)
m = particle_list[0].mass
v = np.linspace(0,10,120)
fv = m*np.exp(-m*v**2/(2*T*k))/(2*np.pi*T*k)*2*np.pi*v
# hist.plot(v,fv, label = "Maxwell–Boltzmann distribution") 
# hist.legend(loc ="upper right")


# from matplotlib.widgets import Slider

# slider_ax = plt.axes([0.1, 0.05, 0.8, 0.05])
# slider = Slider(slider_ax,      # the axes object containing the slider
#                   't',            # the name of the slider parameter
#                   0,          # minimal value of the parameter
#                   tfin,          # maximal value of the parameter
#                   valinit=0,  # initial value of the parameter 
#                   color = '#5c05ff' 
#                  )

# def update(time):
#     i = int(np.rint(time/timestep))
    
#     #ax.set_title('Energy =' + str(Energy[i]))
    
#     # Draw Particles as circles
#     for j in range(particle_number):
#         circle[j].center = particle_list[j].solpos[i][0], particle_list[j].solpos[i][1]
#     hist.clear()    
    
#     # Graph Particles speed histogram
#     vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]
#     hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
#     hist.set_xlabel("Speed")
#     hist.set_ylabel("Frecuency Density")
    
#     # Compute 2d Boltzmann distribution
#     E = total_Energy(particle_list, i)
#     Average_E = E/len(particle_list) 
#     k = 1.38064852e-23
#     T = 2*Average_E/(2*k)
#     m = particle_list[0].mass
#     v = np.linspace(0,10,120)
#     fv = m*np.exp(-m*v**2/(2*T*k))/(2*np.pi*T*k)*2*np.pi*v
#     hist.plot(v,fv, label = "Maxwell–Boltzmann distribution") 
#     hist.legend(loc ="upper right")
    
# slider.on_changed(update)
# plt.show()

Posiciones x

In [23]:
pos_x_list = [particle_list[j].solpos[0][0] for j in range(len(particle_list))]

pos_x_list = np.array(pos_x_list)

pos_x_list.shape
Out[23]:
(300,)
In [24]:
pos_x_hist = np.histogram(pos_x_list, bins=10, range=None)

pos_x_hist
Out[24]:
(array([36, 27, 25, 28, 36, 28, 26, 31, 27, 36]),
 array([ 1.5450174 , 11.38373108, 21.22244476, 31.06115845, 40.89987213,
        50.73858581, 60.57729949, 70.41601317, 80.25472685, 90.09344054,
        99.93215422]))
In [25]:
entropy_x = np.sum( pos_x_hist[0]*np.log(pos_x_hist[0]) )
entropy_x
Out[25]:
1023.2346976366081

Posiciones y

In [26]:
pos_y_list = [particle_list[j].solpos[0][1] for j in range(len(particle_list))]
In [27]:
pos_y_hist = np.histogram(pos_y_list, bins=10, range=None)

pos_y_hist
Out[27]:
(array([37, 29, 24, 28, 31, 31, 30, 26, 33, 31]),
 array([ 1.64801777, 11.45382488, 21.259632  , 31.06543912, 40.87124623,
        50.67705335, 60.48286047, 70.28866758, 80.0944747 , 89.90028182,
        99.70608893]))
In [28]:
entropy_y = np.sum( pos_y_hist[0]*np.log(pos_y_hist[0]) )
entropy_y
Out[28]:
1022.3225510346361
In [52]:
particle_list[0].solpos[0:5]
Out[52]:
[array([11.274827  , 83.31879714]),
 array([11.10109659, 83.23023788]),
 array([10.92736617, 83.14167862]),
 array([10.75363576, 83.05311935]),
 array([10.57990534, 82.96456009])]
In [31]:
pos_y_hist[0]
Out[31]:
array([37, 29, 24, 28, 31, 31, 30, 26, 33, 31])
In [32]:
entropies_y = []

for k in range(len(particle_list[0].solpos)):
    pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
    pos_y_hist = np.histogram(pos_y_list, bins=50, range=(0,200))
    entropy_y = -np.sum( pos_y_hist[0]*np.log(pos_y_hist[0]+1) )
#     entropy_y = -np.sum( np.log(pos_y_hist[0]**pos_y_hist[0]) )
    entropies_y.append(entropy_y)

entropies_y = np.array(entropies_y)
entropies_y[0:5]
#     print(entropy_y)
Out[32]:
array([-776.14128651, -777.32283156, -776.76383835, -776.57980253,
       -776.38504003])

Entropia

Asumiendo probabilidades independientes para $x$, $y$, $p_x$ y $p_y$

$ p = p_x \ p_y \ p_{px} \ p_{py} $

$ S = -p \ln{p} \\ S = -p_x \ln{p_x}-p_y \ln{p_y}-p_{px} \ln{p_{px}}-p_{py} \ln{p_{py}} $

$ S = S_x + S_y + S_{px} + S_{py} $

Entropia de $y$

$ S_y = -p_y \ln{p_y} $

In [55]:
plt.figure(figsize=(12,7))
plt.plot(entropies_y)
Out[55]:
[<matplotlib.lines.Line2D at 0x7f19887c8940>]
In [34]:
entropies_x = []

for k in range(len(particle_list[0].solpos)):
    pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
    pos_x_hist = np.histogram(pos_x_list, bins=50, range=(0,200))
    entropy_x = -np.sum( pos_x_hist[0]*np.log(pos_x_hist[0]+1) )
#     entropy_x = -np.sum( np.log(pos_x_hist[0]**pos_x_hist[0]) )
    entropies_x.append(entropy_x)

entropies_x = np.array(entropies_x)
entropies_x[0:5]
#     print(entropy_y)
Out[34]:
array([-772.5424401 , -772.22401752, -773.32299829, -773.17725408,
       -772.61153761])
In [62]:
# plt.figure(figsize=(12,7))
# plt.plot(entropies_x)
In [61]:
# plt.figure(figsize=(12,7))
# plt.plot(entropies_x, entropies_y)

Entropia posicion

In [58]:
plt.figure(figsize=(12,7))
plt.plot(entropies_x + entropies_y)
Out[58]:
[<matplotlib.lines.Line2D at 0x7f19886913d0>]
In [38]:
vel_x_list = [particle_list[j].solvel[15][0] for j in range(len(particle_list))]
vel_x_list[0:5]
Out[38]:
[-2.6727756063889205,
 -1.867853846575191,
 -1.7728784576290662,
 -2.9653351172118225,
 -2.982865902808966]
In [39]:
vel_x_hist = np.histogram(vel_x_list, bins=10, range=(0,10))
vel_x_hist
Out[39]:
(array([5, 5, 7, 0, 0, 0, 0, 0, 0, 0]),
 array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.]))
In [40]:
entropies_vx = []

for k in range(len(particle_list[0].solpos)):
    vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
    vel_x_hist = np.histogram(vel_x_list, bins=50, range=(-5,5))
    entropy_vx = -np.sum( vel_x_hist[0]*np.log(vel_x_hist[0]+0.01) )
    entropies_vx.append(entropy_vx)
    
#     vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]

#     vel_hist = np.histogram(vel_mod, bins=10, range=None, normed=None, weights=None, density=True)
    
entropies_vx = np.array(entropies_vx)
entropies_vx[0:5]
#     print(entropy_y)
Out[40]:
array([-962.05608072, -958.14838194, -958.14838194, -958.14838194,
       -944.25461567])
In [41]:
entropies_vy = []

for k in range(len(particle_list[0].solpos)):
    vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]
    vel_y_hist = np.histogram(vel_y_list, bins=50, range=(-5,5))
    entropy_vy = -np.sum( vel_y_hist[0]*np.log(vel_y_hist[0]+0.01) )
    entropies_vy.append(entropy_vy)
    
#     vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]

#     vel_hist = np.histogram(vel_mod, bins=10, range=None, normed=None, weights=None, density=True)
    
entropies_vy = np.array(entropies_vy)
entropies_vy[0:5]
#     print(entropy_y)
Out[41]:
array([-937.30635894, -937.30635894, -937.30635894, -932.28208167,
       -932.28208167])
In [42]:
 entropies_vy
Out[42]:
array([-937.30635894, -937.30635894, -937.30635894, ..., -602.93158629,
       -610.86101268, -606.53632224])
In [43]:
np.histogram(vel_y_list, bins=10)
Out[43]:
(array([ 9, 11, 34, 51, 73, 67, 35, 16,  2,  2]),
 array([-5.94319067, -4.68727709, -3.4313635 , -2.17544992, -0.91953634,
         0.33637724,  1.59229083,  2.84820441,  4.10411799,  5.36003157,
         6.61594516]))
In [59]:
Entropia_r = entropies_x + entropies_y #+ entropies_vx + entropies_vy
Entropia_v = entropies_vx + entropies_vy

plt.figure(figsize=(12,7))
plt.plot(Entropia_r)
plt.show()

plt.figure(figsize=(12,7))
plt.plot(Entropia_v)
plt.show()

Entropia = Entropia_r + Entropia_v

Entropia

In [60]:
plt.figure(figsize=(12,7))

plt.plot(Entropia)
plt.show()

Figura para el video

In [ ]:
# for k in range(stepnumber):
# for k in range(0, stepnumber, 50):

k = 1#stepnumber

fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(2,2,1)

# hist = fig.add_subplot(1,2,2)

plt.subplots_adjust(bottom=0.2,left=0.15)

ax.axis('equal')
ax.axis([-1, 30, -1, 30])

ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_xlim([0,boxsize])
ax.set_ylim([0,boxsize])

# Draw Particles as circles
circle = [None]*particle_number
for i in range(particle_number):
    circle[i] = plt.Circle((particle_list[i].solpos[k][0], particle_list[i].solpos[k][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
    ax.add_patch(circle[i])

vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]    
vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]

pos_partic_j = np.array(particle_list[31].solpos)

plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])

pos_partic_j = np.array(particle_list[75].solpos)

plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])

pos_partic_j = np.array(particle_list[42].solpos)

plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])

ax.quiver(pos_x_list, pos_y_list, vel_x_list, vel_y_list, color='orange')

hist = fig.add_subplot(2,2,2)

vel_mod = [particle_list[j].solvel_mag[k] for j in range(len(particle_list))]
hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
hist.set_xlabel("Speed")
hist.set_ylabel("Frecuency Density")

# Compute 2d Boltzmann distribution
E = total_Energy(particle_list, k)
Average_E = E/len(particle_list) 
kb = 1.38064852e-23
T = 2*Average_E/(2*kb)
m = particle_list[0].mass
v = np.linspace(0,10,120)
fv = m*np.exp(-m*v**2/(2*T*kb))/(2*np.pi*T*kb)*2*np.pi*v
hist.plot(v,fv, label = "Maxwell–Boltzmann distribution") 
hist.legend(loc ="upper right")

hist_x = fig.add_subplot(2,2,3)

pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]

hist_x.hist(pos_x_list, bins= 50, range=(0,200), density = True, label = "Simulation Data")

hist_y = fig.add_subplot(2,2,4)

pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]

hist_y.hist(pos_y_list, bins= 50, range=(0,200), density = True, label = "Simulation Data")

# plt.quiver(x_pos, y_pos, final_vector[0], final_vector[1], scale=1e9, units="xy")

Hago las imagenes para el video (esto tarda como 10 minutos)

In [63]:
fig, ax = plt.subplots(figsize=(12, 12))

for k in range(0, stepnumber, 5): #stepnumber
    
    ax.clear()
    plt.clf()
    # plotting line
#     ax.plot(x, np.sinc(x**2) + np.sin(x + 2 * np.pi / duration * t), lw = 3)
#     ax.set_ylim(-1.5, 2.5)

    # returning numpy image
#     return mplfig_to_npimage(fig)


    # hist = fig.add_subplot(1,2,2)

    # plt.subplots_adjust(bottom=0.2,left=0.15)

    
    ax = fig.add_subplot(2,2,1)

    # hist = fig.add_subplot(1,2,2)

    plt.subplots_adjust(bottom=0.2,left=0.15)

    ax.axis('equal')
    ax.axis([-1, 30, -1, 30])

    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)
    ax.set_xlim([0,boxsize])
    ax.set_ylim([0,boxsize])

    # Draw Particles as circles
    circle = [None]*particle_number
    for i in range(particle_number):
        circle[i] = plt.Circle((particle_list[i].solpos[k][0], particle_list[i].solpos[k][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
        ax.add_patch(circle[i])

    vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]    
    vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
    pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
    pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]

    pos_partic_j = np.array(particle_list[31].solpos)

    plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])

    pos_partic_j = np.array(particle_list[75].solpos)

    plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])

    pos_partic_j = np.array(particle_list[42].solpos)

    plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])

    ax.quiver(pos_x_list, pos_y_list, vel_x_list, vel_y_list, color='orange')

    hist = fig.add_subplot(2,2,2)

    vel_mod = [particle_list[j].solvel_mag[k] for j in range(len(particle_list))]
    hist.hist(vel_mod, bins= 12, density = True, label = "Simulation Data")
    hist.set_xlabel("Speed")
    hist.set_ylabel("Frecuency Density")

    # Compute 2d Boltzmann distribution
    E = total_Energy(particle_list, k)
    Average_E = E/len(particle_list) 
    kb = 1.38064852e-23
    T = 2*Average_E/(2*kb)
    m = particle_list[0].mass
    v = np.linspace(0,10,120)
    fv = m*np.exp(-m*v**2/(2*T*kb))/(2*np.pi*T*kb)*2*np.pi*v
    hist.plot(v,fv, label = "Maxwell–Boltzmann distribution") 
    hist.legend(loc ="upper right")
    
    plt.ylim([0, 0.4])

    
    hist_x = fig.add_subplot(2,2,3)

    pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]

    hist_x.hist(pos_x_list, bins= 15, range=(0,200), density = True, label = "Simulation Data")
    
    plt.title('Histograma de x')
    
    plt.ylim([0, 0.02])

#     hist_y = fig.add_subplot(2,2,4)

#     pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]

#     hist_y.hist(pos_y_list, bins= 20, range=(0,200), density = True, label = "Simulation Data")

#     plt.ylim([0, 0.02])

    entrop_plot = fig.add_subplot(2,2,4)

    entrop_plot.plot(Entropia[:k])
    
    plt.title('Entropía')
    
    plt.ylim([Entropia.min(), Entropia.max()])
    
    plt.xlim([0, stepnumber])
    
    plt.xlabel('Pasos de tiempo')
    
    plt.savefig(f'gas_t_{k}.png', dpi=100)
In [48]:
!pip install moviepy

Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: moviepy in /home/hugo/.local/lib/python3.8/site-packages (1.0.3)
Requirement already satisfied: tqdm<5.0,>=4.11.2 in /home/hugo/.local/lib/python3.8/site-packages (from moviepy) (4.64.0)
Requirement already satisfied: imageio<3.0,>=2.5 in /home/hugo/.local/lib/python3.8/site-packages (from moviepy) (2.16.1)
Requirement already satisfied: requests<3.0,>=2.8.1 in /usr/lib/python3/dist-packages (from moviepy) (2.22.0)
Requirement already satisfied: decorator<5.0,>=4.0.2 in /usr/lib/python3/dist-packages (from moviepy) (4.4.2)
Requirement already satisfied: numpy in /home/hugo/.local/lib/python3.8/site-packages (from moviepy) (1.20.3)
Requirement already satisfied: proglog<=1.0.0 in /home/hugo/.local/lib/python3.8/site-packages (from moviepy) (0.1.9)
Requirement already satisfied: imageio-ffmpeg>=0.2.0 in /home/hugo/.local/lib/python3.8/site-packages (from moviepy) (0.4.5)
Requirement already satisfied: pillow>=8.3.2 in /home/hugo/.local/lib/python3.8/site-packages (from imageio<3.0,>=2.5->moviepy) (9.1.0)
WARNING: You are using pip version 21.3.1; however, version 22.0.4 is available.
You should consider upgrading via the '/usr/bin/python3 -m pip install --upgrade pip' command.

Junto las imagenes en el video

In [64]:
from moviepy.editor import ImageSequenceClip

clip = ImageSequenceClip([f'gas_t_{k}.png' for k in range(0, stepnumber, 5)], fps=20) #stepnumber

clip.ipython_display(fps = 20, loop = True, width=720)

Moviepy - Building video __temp__.mp4.
Moviepy - Writing video __temp__.mp4


                                                              

Moviepy - Done !
Moviepy - video ready __temp__.mp4
Out[64]: